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We provide detailed analysis of the complex energy eigenvalue spectrum for a two-channel quan- 
tum wire with an attached adatom impurity. The study is based on our previous work [Phys. Rev. 
Lett. 99, 210404 (2007)], in which we presented the quasi-bound states in continuum (or QBIC 
states) . These are resonant states with very long lifetimes that form as a result of two overlapping 
continuous energy bands one of which, at least, has a divergent van Hove singularity at the band 
edge. We provide analysis of the full energy spectrum for all solutions, including the QBIC states, 
and obtain an expansion for the complex eigenvalue of the QBIC state. We show that it has a small 
decay rate of the order g 6 , where g is the coupling constant. As a result of this expansion, we find 
that this state is a non-analytic effect resulting from the van Hove singularity; it cannot be predicted 
from the ordinary perturbation analysis that relies on Fermi's golden rule. We will also numerically 
demonstrate the time evolution of the QBIC state using the effective potential method in order to 
show the stability of the QBIC wave function in comparison with that of the other eigenstates. 

PACS numbers: 



I. INTRODUCTION 

In a previous Letter Q we introduced the quasi-bound 
states in continuum (QBIC) as resonant states which oc- 
cur in certain systems with overlapping continuous en- 
ergy bands. Consider the case where one of these en- 
ergy bands has a divergent van Hove singularity in the 
density of states at one of the overlapping band edges. 
Then a discrete excited state coupled to these energy 
bands gives rise to a metastable resonant state with an 
extended lifetime. This effect cannot be predicted using 
Fermi's golden rule as it breaks down in the vicinity of 
the singularity. 

In our Letter, we demonstrated the existence of the 
QBIC states in the context of a two-channel quantum 
wire coupled to a single adatom impurity. This built on 
previous work by S. Tanaka, S. G., and T. P. on a sin- 
gle channel wire coupled to an adatom [2J , in which they 
demonstrated various non-analytic effects that resulted 
from the presence of the divergent van Hove singularity 
in the electron density of states (DOS) [3| at the edge of 
the conduction band; note that these are characteristic 
effects of the van Hove singularity in a one-dimensional 
system. One of these effects was a bound state that lies 
just outside of either edge of the conduction band, no 
matter how deeply the discrete adatom energy is em- 
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bedded in the band @, H H H|. We will refer to this 
state here as a persistent bound state as it co-exists with 
the unstable decay states. In the present paper, we will 
show that when the adatom is coupled to a second en- 
ergy band that overlaps the first, this persistent bound 
state is slightly destabilized due to the fact that it lies in 
the continuum of the second energy band. We will ex- 
plicitly demonstrate that it is this persistent bound state 
(when de-stabilized) that forms the QBIC eigenstate for 
the two-channel system; for this purpose, we compare 
term-by-term the analytic expansions of the bound state 
energy eigenvalue in the single channel system and that 
of the QBIC energy eigenvalue in the two-channel sys- 
tem. We will also show that the decay rate for the QBIC 
state is (to the lowest order) proportional to g 6 , where 
g is the coupling constant between the adatom and the 
site to which it is attached. We assume that g is small in 
this paper. 

We have envisioned the QBIC state as a generalization 
of the bound state in continuum (BIC) originally pro- 
posed by von Neumann and Wigner in 1929 Since 
their initial proposal, a good deal of theoretical study 
has been devoted to this phenomenon [j| 0, E3, El, 0, 
EE EE EE EH , including a recent paper by S. Longhi, in 
which the author demonstrates the presence of BIC states 
in a single level semi-infinite Fano- Anderson model [l7| 
and another article by S. Tanaka, S. G., G. Ordonez, 
and T. P., in which the authors demonstrate the pres- 
ence of BIC states in a single channel quantum wire 
coupled with multiple adatom impurities [y, [18| . There 
has also been experimental confirmation of the BIC phe- 
nomenon (l9l |2Q |. However, since it is a zero measure 
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effect (meaning that it only occurs at discrete points in 
parameter space), it is generally considered difficult to 
detect. 

As discussed in our previous Letter, the QBIC state 
actually has a decay rate (due to the imaginary part of 
the complex eigenenergy) and hence does not technically 
lie in the continuous energy spectrum of the conduction 
band. However, the decay rate for this state is on the or- 
der of <7 6 , much smaller than the ordinary decay rate that 
is predicted to be of the order g 2 by Fermi's golden rule. 
Hence the QBIC state will behave as a bound state (with 
real part of the complex energy inside the continuous en- 
ergy spectrum) even on relatively large time scales. We 
will also show below that the g 6 power in the decay rate 
is a direct result of the interaction between the divergent 
van Hove singularity at one of the band edges and the 
continuum of the other band. 

While the BIC states occur only at discrete values 
in parameter space, the QBIC states occur over a wide 
range of parameter space. Hence, it may be much easier 
to experimentally verify the QBIC effect. It may also 
be easier to verify the QBIC effect, considering that the 
effect will likely occur in other physical systems; for in- 
stance, in another Letter d by C.-O. Ting, T. P. and 
S. G., the authors explored the effects of the divergent 
van Hove singularity in the photon density of states at the 
cutoff frequency in the interaction between an excited os- 
cillator (diatomic molecule, for instance) and the lowest 
Transverse Electric (TE) mode in a rectangular waveg- 
uide. The QBIC effect will also occur in this waveguide 
system when the oscillator interacts with the second low- 
est TE mode, which has a cutoff frequency embedded in 
the continuum of the lowest TE mode. 

For the present case, we consider the system shown in 
Fig. [3a), which is composed of two tight-binding chains 
and an adatom or quantum dot. The two tight-binding 
chains (labeled y = 1, 2 in Fig. [IJa)) both have internal 
hopping parameter —t^/2 (internal sites of both chains 
are labeled by integer x, where \x\ < m with N = 2m + 1 
and N (^> 1) is the number of sites in either chain y = 1 
or 2). The two chains are then coupled together site-by- 
site with hopping parameter — t' h , creating a ladder shape. 
The dot (labeled d) is then coupled to the x = site of 
the y = 1 chain. Hence, we can write the Hamiltonian 
for our system as 

-^(M)(x,2| + |x,2)(x,l|) 
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FIG. 1: (a) An adatom (quantum dot) attached to a ladder, 
(b) After partial diagonalization in the y direction, the system 
is composed of the dot coupled to two independent channels. 



vice versa). The third term then represents hopping be- 
tween the ad-atom and the (0, 1) site of the ladder and 
finally the fourth term gives the unshifted energy of the 
ad-atom. 

In order to diagonalize the second term of the Hamil- 
tonian ([1]), we introduce the basis 



\x,+)\ = \x,l) 
\x,-) - U \\x,2) 



where 



U=±(] 1 
M 1 1 



= u~ 



(2) 



(3) 



Using the new basis \x, a = ± 
divided into the a — + and — 
term) as 



, the Hamiltonian can be 
chains (with the adatom 



H 



E 



"T ^2(\ x + 1 >°')< a; > (T l + \x,a)(x+ l,a\ 



-25 d |d)<d|. 



(4) 



- 5 (|d)(0,l| + |0,l)(d|) 
-£d|d)(d|, 



(1) 



in which Ed denotes the energy of the dot. In accor- 
dance with our designations in Fig. [T] the first term here 
represents internal hopping along either of the chains (in 
the x direction) while the second term describes hop- 
ping from one chain to the other (y = 1 to y = 2 and 



We have now obtained the Hamiltonian in the form of 
Fig. [Tf b) , in which the two a = +, — chains represent 
two independent channels for charge transfer. Note that 
we can also interpret the label a as electron spin and t' h 
as a magnetic field. 

In the next section, we will outline two approaches to 
obtaining the full diagonalization of the above Hamil- 
tonian. In the first approach, we will introduce the 
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wave vector representation to obtain a bi-linear form 
of the Hamiltonian, from which we can obtain the en- 
ergy eigenvalues of the system using the analysis due to 
Friedrichs [2lJ. In the second approach we will rely on 
the recently presented method of outgoing waves [22| . 

In Sec. Ill we will present the full energy eigenvalue 
spectrum for the case in which the energy bands associ- 
ated with the two channels of the quantum wire overlap; 
we show the energy shift and decay rate of the QBIC 
states. In particular, we will show that the decay rate of 
the QBIC states is proportional to g 6 and that this effect 
is a direct result of the interaction between the van Hove 
singularity at the edge of one energy band and the con- 
tinuum of the other energy band. We will also examine 
the wave function and time evolution for the QBIC state. 

In Sec. IV we will examine the energy spectrum in two 
special cases. In the first case the two channels become 
decoupled (that is, t' h — 0). The energy spectrum then 
reduces to that of a single channel quantum wire coupled 
with an adatom [2|. In the second case the lower edge 
of the upper band coincides with the upper edge of the 
lower band (that is, t' h = t^). We will then examine the 
energy spectrum and discover a modification of the QBIC 
states resulting from the two overlapping singularities. 

Finally, in Sec. V we will briefly outline our results and 
make our final conclusions. We will also discuss the gen- 
eralization to an n-channel model briefly. In Appendix, 
we summarize a numerical method of following the time 
evolution of decaying resonant states. 
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FIG. 2: The two continuous dispersions relations (0 which 
form two conduction bands (channels) E± in the wire system. 
Here we graph the overlapping case in which < t' h < t^. 

vector representation with wave vectors K± in the two 
respective channels a = ± by 



\K+) = 



TTxfE' 



t iK±x 



x,±), 



(5) 



where K± = n±Ak with Ak = 2tt/N and integers n±. 
This allows us to write the Hamiltonian ([4]) as a variation 
of the bi-linear Friedrichs-Fano model 



n 



EE 

<r=± K a 

+£ d |d)(d| 



E a \K a )(K a \ 



\K*)(d\) 
(0) 



II. DISPERSION RELATION AND 
DIAGONALIZATION OF THE HAMILTONIAN 

In this section we will outline two methods for diago- 
nalizing the Hamiltonian J4|. In the first we will write 
the Hamiltonian in a bi-linear form, from which we can 
immediately obtain the eigenvalues following the method 
due to Friedrichs f2l[. These eigenvalues will be ob- 
tained in the form of discrete solutions to a dispersion 
relation that is equivalent to a 12th order polynomial. 
The discrete solutions of this polynomial give the diag- 
onalized energy shifts and decay rates for an electron in 
the adatom. 

In the second approach we will rely on the method 
of outgoing waves presented recently by N. H., H. N., 
K. Sasada, and T. P. [22|. This method will also yield 
the dispersion relation, but is better suited for performing 
numerical simulations. 



A. Dispersion relation from the Friedrichs solution 
for bi-linear Hamiltonian 

Because we are interested in the case N ^ 1, we may 
impose boundary conditions for our convenience of the 
calculations. Imposing here the usual periodic boundary 
conditions in the x direction, we may introduce the wave 



The energies E a in the two channels are determined by 
their respective wave numbers K a according to 



E± = — < h cos K± =F t' h . 



(7) 



We refer to the above equations as the continuous dis- 
persion equations for the system; this is because in the 
continuous limit N — > 00 they describe the allowed ener- 
gies for the two continua of K a states. In Fig. [5] we graph 
these two energy bands for the case t' h < <h, in which 
they will overlap. (By contrast, below we will write a 
discrete dispersion equation that describes how the dis- 
crete energy E<± is modified by the interaction.) Both 
of the continuous channels described in Eq. ([7]) have an 
associated density of states (DOS) function. The DOS 
functions for the two channels are given by 



1 



P±{E) = - 



1 



n^/ti~(E±t' h ) 



(8) 



Note the presence of two van Hove singularities in either 
channel. These singularities are located at ±ih + t' h for 
the "+" channel and ±th — t' h for the "— " channel. 

Since the Hamiltonian ([5]) is in a bi-linear form, in 
principle we may now diagonalize to explicitly solve the 
problem according to the method given by Friedrichs [2l| . 
However, we will leave the issue of obtaining the ex- 
plicit solutions to the following subsection, in which the 
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method is more suited to conducting numerical simula- 
tions for the time evolution of the system. Instead, we 
will simply follow the standard method to obtain the dis- 
crete dispersion relation for an electron inside the two- 
channel quantum wire. 

The discrete dispersion relation is given by 77(2) = 2 — 
E d — = 0, where the self-energy for an electron 
in the two-channel wire is determined by 



E(z) 



2N ^-f ^ 2 - E k 



9 

2N 



E 



z + t' h — th cos fc + 



y I 

z - t' h - t b COS 



4w 



<r_ 

2 



— t[, — Ui cos k. 

1 

z + ti — th cos fc+ 

1 

z — t' h — th cos k- 

1 

+ 



(9) 



Thus we find the dispersion equation 



z- Ea- 



st 
2 



1 



+ 



1 



= 



~ (io) 

as reported previously [l|, la]. This dispersion equation 
describes the behavior of an electron initially trapped in 
the adatom and can be written as a 12th order poly- 
nomial equation in z. Hence, we will also refer to this 
equivalent equation as the dispersion polynomial. The 
twelve discrete solutions to this equation give the allowed 
bound states (purely real solutions) and resonant states 
(complex solutions) in the diagonalized system. As we 
will discuss below, these solutions can be viewed as liv- 
ing in a complex energy surface, parameterized by the 
original impurity energy E d . In the present case of the 
two-channel wire, this energy surface will be composed 
of four Riemann sheets. 

Once we have obtained the twelve discrete solutions 
2 = E to the dispersion polynomial, then Eq. ((7| implies 
that each eigenvalue E can be assigned two K± values. 
With the values K± in hand, we will be able to write the 
wave function for each solution, making use of Eq. (fl~2]) 
given below. Practically, this means that any electron 
state in the wire can be fully described by the three values 
E, K + and K— 



B. Solutions of Hamiltonian by the method of two 
outgoing waves 

We will now solve the Schrodinger equation 

H\il>)=E\1>), (11) 

for the resonant states of the Hamiltonian given in 
Eq. ([!]). As has been previously shown 22, 29], the reso- 
nant cigenfunction of the tight binding model on a chain 
with an adatom can be written in the form 



or 



(12) 



(13) 

Using the resonant eigenfunction (|12[) in the Schrodinger 
equation (fTTj) for the case x ^ 0, we obtain 

Eip a (x) = H ipo-(x) 

= ~{^{x + l)+i) a {x-l)}-at^ a {x) 

= (-t h cos K a - at' h )^ a (x). (14) 

This again yields the dispersion relations E — 
—thCOsK± =p t' h inside the two channels of the wire as 
given in Eq. ([7]) above. 

To solve the eigenequation 

ip d = (d\tp) and ip(x, y) = (x, y\ip), (15) 

for x = and x = d, respectively, we return to the orig- 
inal y space. Using the original bases \x,y) and Eq. (J2J, 
the resonant eigenfunctions i/j(x, y) in the original y space 
have the form 



ip y (x) 



i>(x,l) 
= UMx) 

-;W™(D + > jr - w (- 1 > 

(16) 

With the Hamiltonian (p} the Schrodinger equa- 
tions (fTTj) for x — and x — d become 



■| {V>(-1, 1) + V(U 1)} - t' h ^(0, 2) + g^ d 



i?V(o,i), 

^(0,2), 



-yW(-l, 2) +V(1,2)} -^(0,1) 

^(0,1) + E d ^ d = E^j d . 

Substituting the resonant wave functions ip in the site 
representation from Eq. (|16p while making use of the con- 
tinuous dispersion relations Eq. (O, we obtain 

ithA+ sin i4T+ + ithA- sin K- - V2gip d = 0, 
it h A + sin K + - it h A_ sin K _ = 0, (17) 
g {A+ + A-) + s[2 (E d - E) V> d = 0, 
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which can be written in matrix form as 

ith sin K+ ith sin K- — \/2<? \ { A+\ 
it h sinK + -zi h sinK_ I A_ J = 0. 

g g V2(E d ~E)J \Vd / 

(18) 

In order to have non-trivial solutions to Eq. (| 18[> . the 
determinant of the coefficient matrix above must be zero. 
Hence we obtain the following condition on A + , A_ and 
ipd- 



E - E c \ 



1 



2iih sin K + 2ith sin K- 



(19) 



Making use of the continuous dispersion equations (|7|). 
we can see that the above condition is equivalent to the 
discrete dispersion equation (fTO)) related to the interac- 
tion between the adatom and the two a channels. 

According to the previous work 0, H3] , the dispersion 
equation for the single-chain model with an adatom is 
given by 



E 



cha 



Ea 



2ith sin K r 



(20) 



Hence we note that the present dispersion equation (|T9[) 
for the ladder model corresponds to the sum of two single- 
chain dispersion equations. 

In the case of the single chain model 0, 0, [22|, the 
complex energy spectrum could be evaluated in terms of 
a complex plane consisting of two Riemann sheets. In 
that case there was only one wave number -KT c h a in corre- 
sponding to the single channel available to an electron. 
One can then easily classify whether a resonant state lies 
in the first or second Riemann sheet according to the 
sign of the imaginary component of this wave number. 
This determination is beneficial in obtaining a detailed 
understanding of the survival probability for the excited 
state, in which a contour deformation must be performed 
in the complex energy plane. For instance, the position 
of the poles may also influence the strength of the non- 
Markovian effect due to the so-called branch point effect 
(this will be the subject of a future publication). 

In the present case of the two-channel model, there are 
two wave numbers K± resulting from the two channels 
available to the electron. The imaginary part of these 
two wave numbers together provides four possible sign 
combinations and hence the complex energy plane is now 
composed of four Riemann sheets; see Fig. [3] for an ex- 
ample in the case < t' h < t^. We can also see that the 
Riemann surface must be four-sheeted by considering the 
dispersion equation (|10p , in which each of the two roots 
may take either a positive or negative sign, again result- 
ing in four combinations (although one must be careful 
here as the sign combinations in this approach change 
for different portions of the same sheet). In the general 
case of an rt-channel quantum wire, the complex energy 
surface will be composed of 2n Riemann sheets. We will 
comment further on the generalization to an arbitrary 
number of channels later in this paper. 




FIG. 3: Four sheeted Riemann energy surface for the two 
channel model in the case < t' h < th- Solutions of 
the discrete dispersion equation with the sign combination 
(sgn(Imif+), sgn(Imif-)) = (+, +) lie in the first Riemann 
sheet (solid black line). Solutions with the combination ( — , +) 
he in Sheet II (short-dashed green line), those with the com- 
bination (+, — ) line in Sheet III (long-dashed red line), and 
those with the combination (— ,— ) lie in Sheet IV (chained 
blue line). The curved lines in the center of the diagram rep- 
resent the two branch cuts where the sheets intersect along 
their respective real axes. For instance, on the left side of 
the diagram, if one starts from the positive imaginary half 
of Sheet I and crosses the real axis between — tj, — t' h and 
— th + t' h (represented by the curved overlapping solid black 
and short-dashed green lines) then one will emerge on the 
negative imaginary half of Sheet II. 



For the purpose of assignment of each solution to the 
correct Riemann sheet, it is more convenient to mod- 
ify Eq. (fT9|) with the help of the continuous dispersion 
equation (O and solve the following set of simultaneous 
equations with respect to K± than to solve the discrete 
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I ! LI tmt* Hp hUhl*4f< 
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FIG. 4: Band structure for the case < t' h < t h . See Eq. J7J) 
for the definition of _E± . 

dispersion equation (|10p with respect to z = E: 
—th cos K+ — t' h = — th cos K _ + ij, 

^ d+5 (2^ sin if + + 2it h sini<:_) ' ^ 

The signs of Im K ± of each solution give the correct Rie- 
mann sheet immediately. 

Note that the first Riemann sheet is assigned in a nat- 
ural way as the energy eigenvalues of all of the solutions 
in this sheet must be real, each corresponding to a bound 
state in the energy spectrum analysis. (They are on the 
positive imaginary axes of the wave-number spaces K±.) 
This is because the Hamiltonian must behave in a man- 
ner equivalent to a Hermitian operator in this first sheet 
with all eigenvalues being purely real. It is only when the 
Hamiltonian is extended into the rigged Hilbert space [23[ 
that complex solutions in the other sheets may be con- 
sidered on an equal footing with the stable solutions in 
the first sheet, and that they can be interpreted as com- 
plex eigenvalues of the Hamiltonian in the rigged Hilbert 
space |2H [25|. 

III. ENERGY SPECTRUM ANALYSIS FOR 
< t' h < t h CASE AND QBIC STATES 

In this section we analyze the eigenenergy spectrum 
for the case < t h < th, in which the two conduction 
bands overlap. In this overlapping case the edge of one 
band is embedded in the continuum of the other and vice- 
versa. This results in two outer band edges and two inner 
(embedded) band edges (see Fig. 2]). 

When the energy of the impurity lies inside the overlap- 
ping region between the two inner band edges (— th+t h < 
Ed < th — t'h), we will find that there exist persistent sta- 
ble solutions that lie just outside of the two outer band 
edges, as mentioned in Introduction 0, 0, H, @. As we 
will see, these stable states are a direct result of the van 
Hove singularity in the density of states at the band edge. 
For the inner band edges, however, there will now be two 
competing effects: the first being the stabilizing effect of 
the singularity from the embedded band edge and the 
second being the de-stabilizing effect of the continuum in 
which it is embedded (i.e., the second conduction band). 
We will see that this will result in a slightly de-stabilized 
state embedded in the continuum (the QBIC state). 



A. Energy spectrum analysis from the dispersion 
polynomial 

We will now analyze in detail the twelve solutions to 
the discrete dispersion equation (fTUj) . Before presenting 
the full complex energy spectrum as a function of Ed , we 
will first consider the eigenenergies at two specific val- 
ues of Ed in Tables Q] and [III m order to illustrate our 
earlier point regarding the placement of the solutions in 
the complex energy surface. In these tables, we present 
the values obtained for each of the twelve solutions by 
solving the dispersion equation (fTU]) for the eigenener- 
gies E and the continuous dispersion equations ([7]) for 
the wave number pairs K± at the two values Ed = 0.3ih 
and Ed = —1.0th, respectively. For the other parameters 
of the system, in both tables we have used the values 
t' h = 0.345th and g = O.lih- In these tables (and through- 
out this paper) we measure energy in units th = 1. In 
addition, in Fig. [5] we have indicated the relative position 
of each individual pole in the complex K + , K- and E 
planes, respectively, each given for the value Ed = 0.3th 
from Table |H 

As mentioned previously, the placement of each solu- 
tion in the complex energy plane can be determined in 
a straight-forward manner by the sign of the imaginary 
parts of the two complex wave vectors K±, which are 
given as the solutions of the simultaneous equations (|2ip ; 
see Fig. [6] We have designated those with positive imag- 
inary K + component and positive imaginary K- compo- 
nent ((+,+), respectively) as lying in Riemann Sheet I; 
likewise we have designated (— , +) as Sheet II, (+, — ) as 
Sheet III and (— , — ) as Sheet IV. The resulting eigenen- 
ergy E — — th cos K± =F t' h of each solution is shown in 
Fig. In Tables Q] and [ft] and Figs. © and we have 
labeled each solution by a letter P, Q, R or S according 
to the Riemann sheet (I, II, III, or IV, respectively) in 
which that solution lies. 

As these solutions are roots of a polynomial with real 
coefficients, each complex decay solution (with negative 
imaginary part) is accompanied by a complex conjugate 
"growth" solution with positive imaginary part. How- 
ever, when we calculate the survival probability for the 
excited impurity state, only the decay solutions will con- 
tribute a pole in the complex contour integration. Hence, 
our focus will be the decay solutions indicated by shades 
in Tables U and HD 

Note that Sheet I contains only purely real solutions. 
This choice is necessary as our Hamiltonian should be de- 
fined to reduce to a Hermitian operator (with real eigen- 
values) in the first sheet in the finite case (closed system) 
or when the coupling vanishes. Also notice that there is 
one conjugate pair of solutions (RQ4 and RQ5) labeled 
by two sheets. This is due to the fact that these solu- 
tions lie in Sheet III for negative values of Ed then cross 
into Sheet II (through the double branch cut on the real 
axis of the complex energy plane) for positive values of 
Ed- At the same time, the imaginary part of the energy 
eigenvalue changes its sign; the solution RQ5 is the decay 
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TABLE I: The twelve discrete eigenvalues for t' h = 0.345th, 9 = 0.1th, and — 0.3th- The decay solutions are indicated by 
shading. 



state 


E/t h 


K+ 


K- 


Riemann Sheet 


Pl(0.3) 
P2(0.3) 


1 34501 1 59 
-1.34500463 


3 14159965 4-4 1 11593956 
+i 0.00304629 


3 14159965 4-4 00480148 
+i 1.11592751 


I 
I 


Ql(0.3) 


1.34501136 
-0.655(11371) -/ 1.5D93 / Id" 7 
i • i i'.-.i i 170 i |.5in:', 'IIP 7 
0.29998854 -jO.iMjl 53774 


3.14159265 -i 1.11593245 
1 25558888 , 1.5875 III" 7 


3.14159265 +i 0.00476787 

0. 28*2 f 0.O0-523531 

2*82 523531 

[ 52576970 i. 0.00 1 53930 


II 


•' 




R(.)5(0.3) 


0.2999s 5 1 - 015 !77 1 


-2.27180290 -/ 0.0020122.1 


1.52570970 -/ 0.00153930 
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Rl(0.3) 


-1.34500459 
0.65509906 i 2.9331 Ml " 


+i 0.00303273 
3.1 113.8429 +; 0.014(17702 


-i 1.11592748 
1.88609355 -; 3.0852 VICT" 


in 


- H:!(0 ; ' 


0.05509900 2.9331 " 10 " 


!.1113f 12!) - 7702 


-1 SSI 09 , -/ !.0852 K) ' 


mi 












S2H>.3i 1 


0.29991927 51171, 


-2.27161773 -i 0.01510419 


-1.525703 !3 -/ 0.0 1 155025 





TABLE II: The twelve discrete eigenvalues for t' h — 0.345th, 9 = 0.1th, and — —1.0th- The decay solutions are indicated by 
shading. 



state 


E/t h 


K+ 


K- 


Riemann Sheet 


Pl(-l) 
P2(-l) 


1.34500228 
- 1.34510721 


3.14159265 +i 1.11592578 
+i 0.01464344 


3.14159265 +i 0.00213553 
+i 1.11600280 


I 
I 


Ql(-l) 

Q2I 11 

Q3(-l) 


1.34500226 

I.IM.54-5.-.7G / i Mil 1659855 

105 15070 +/ 0.00059855 


3.14159265 -i 1.11592577 


3.14159265 +i 0.00212886 
- 0.0072i»999 +; 0.81454288 


II 


1910336 -* 0.00S78757 


0.00726999 +i 0.81454288 


n 


r;i, i, 

R3(-l) 
RQ4(-1) 

BQTn-is 


- 1.34510275 

0.65500456 +i 2.9002 xl0~ 8 

- 0.65510500 +i 3.2049 xlO" 6 

- dammit -t lumtxio " 


4-i 0.01433550 


-i 1.11599952 
I..xs599ll5 ■ 3.0505 1" ' 


in 


3.1 115S305 -/' 1.1.(10302110 
1.25549284 +i 3.3711 xl0~ 6 
- j.2.>5«riMi -i 3.3711x10 " 


* - 1 ss 509 115 -. 3.0505 x 10 - s 
- 0.00022112 -i 0.01449328 
a»HL'21 12 -i Wii Him* 


in 
in 


SLl-L) 


- o.«h:u«)7 -i tummm;m 


UMU12L7 -> 0.0(18726: 57 






S2i - 1 l 


— 0.99434007 +i 0.00663666 


- so 1 1 12 17 -/ » 72037 


- 0.0071.1838 -/ 0.802.1.83 10 


n 



solution for E<$ < 0, but the solution RQ4 is the one for 
Ed > 0. These solutions disappear for the case % = t' h 
as will be discussed in Sec. IV. 

We now analyze the detailed energy spectrum for our 
ladder model in the present < t' h < % case. In Fig. [8] 
we present the real part of the twelve solutions of the 
dispersion equation (|10[) as a function of the impurity 
energy We have also plotted the line Rei? = Ed 
that represents an unshifted energy for the adatom (the 
impurity energy if there was no interaction). Hence, the 
deviation of each solution from this line represents the 
energy shift due to the interaction with the two-channel 
wire. In these and the following figures we will always 
use = 1 as the unit of energy. 

The behavior of the two purely real solutions P2 and 
Rl (Fig. 0(e)) are consistent with the energy spectrum 
previously pointed out 0, H, § for the persistent stable 
states mentioned above. For values of Ed -C — fh — t% 
far below the lowest band edge, these two solutions are 
shifted downwards slightly from the line Rei? = Ed- 
(The shift for both solutions can be shown to be pro- 
portional to g 2 , though P2 always has the slightly larger 
shift.) For values Ed ^ —t^ — tL anywhere above the low- 
est band edge we find these solutions are shifted down- 



wards instead from the lowermost band edge Rei5 = 
— th — t' h , consistent with the previously reported behav- 
ior for the persistent stable state. Below, we will show 
that in this case the shift is proportional to g 4 . 

For values Ed -C — <h — t' h we find that the solutions 
SI and S2 (lower left-hand corner of Fig. [5Jd)) are also 
purely real. These states may be called anti-bound states 
or virtual states [26( . As we increase the value of Ed such 
that Ed < — th — t' h we find that these solutions merge 
abruptly to form a complex conjugate pair. This is simi- 
lar to the behavior of the complex solutions in the single 
channel model § . The imaginary part of these solutions 
can be seen in Fig. [9] (graphed with all the complex solu- 
tions), and in detail in Fig. [TTJJ Notice in these figures 
that the decay rate is amplified for SI and S2 in the vicin- 
ity of each of the four band edges. This amplification is 
a result of the breakdown of Fermi's golden rule in the 
vicinity of the van Hove singularity at each of the band 
edges. Ordinarily, the golden rule predicts that to lowest 
order the decay rate will be proportional to square of the 
coupling constant g 2 ; however, the van Hove singularity 
in the context of a one-dimensional system results in a 
decay rate with a non-analytic dependence on the cou- 
pling constant, such that to lowest order the decay rate 
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FIG. 5: The relative positions of the poles in the complex K+, K- and E planes for the value Ed = 0.3, as well as t' h = 0.345 
and g — 0.1, as given in Table|T] The two vertical dotted lines in the figures (a) and (b) indicate the edges of the Brillouin zone 
KeK± — ±7r, while the four vertical dotted lines in the figure (c) represent the van Hove singularities. 



is proportional to g 4 ^ 3 , as previously reported @, Q. 



B. Van Hove singularity and the origin of the 
quasi-bound states in continuum 

We will now consider the detailed behavior and the ori- 
gin of the QBIC effect. The decaying solutions Q2, R2 
and RQ4 (or RQ5) (and their conjugate "growth" part- 
ners) each display this effect in different regions of the 
energy spectrum. Here we will focus on the solution Q2 
as our example in order to demonstrate the properties 
of the QBIC states. Taking symmetry into account (see 
Fig.[9j), the solution R2 behaves in a manner almost anal- 
ogous to Q2. Meanwhile a detailed analysis of the inte- 
gration contour for the survival probability of the excited 
state reveals that the solutions RQ4 and RQ5 may be of 
less significance in terms of the QBIC effect as it does not 
contribute a pole (exponential decay) in the most natu- 
ral integration contour. However, it is possible that this 
solution will play a more significant role in terms of the 
non-Markovian decay. 

Focusing on the solution Q2, we see in Figs. [8] and [9] 
that the solution is complex for values of the impurity 
energy just below the lower inner band edge tj, — t' h . 
However, in Fig. [5] we see that as we increase the value of 
Ed, the real part of Q2 approaches the inner band edge 
ih — £jj in a manner similar to the persistent stable state 
P2 discussed above. Meanwhile, the imaginary compo- 
nent associated with this solution does not vanish near 
the inner band edge. Instead, the decay rate lies near 
the value zero as can be seen in Fig. [9] and in greater de- 
tail in Fig. QT] This is the QBIC state introduced in our 
previous Letter As can be seen from the figures, the 
real part of the energy is embedded in the lower energy 
band similar to the bound states in continuum proposed 
by von Neumann and Wigner, while the decay rate is 
non-zero but remarkably small. 

We will now show that the QBIC effect is a direct result 
of two competing effects resulting from the embedding of 



the van Hove singularity at the edge of one conduction 
band in the continuum of the other band. The first effect 
is the the tendency of the embedded singularity to cre- 
ate a persistent stable state; in other words, there would 
be an ordinary persistent stable state if it were not for 
the second energy band. The second effect is the ten- 
dency of the embedding conduction band to destabilize 
an otherwise stable state. In order to make this point 
explicit, we will obtain an analytic expansion for the en- 
ergy eigenvalue of the persistent stable state in a single 
channel model (see Fig. [T2")) and compare this term-by- 
term with a similar expansion for the QBIC state Q2 
in the present case (see Fig. IT3l) . This follows from our 
discussion in the previous Letter [l[ . 



1. Analytic approximation for the eigenenergy for the 
persistent stable state in single channel model 

We may write the Hamiltonian 7i_ for a single channel 
quantum wire [27| with energy shifted by t' h as 



H- 



~2 ^ 



(\x + l)(x\ + \x){x + l\)+t' h 



+^(|d}(0,l| + |0,l)(d|)+£ d |d)(d|. (22) 

We have chosen the energy offset t' h here such that the 
single energy band for this Hamiltonian mimics that of 
the upper energy band (— t^cosK- + t' h ) in the two- 
channel model. This single channel will also have the 
same density of states function p_ from Eq. JHJ). Then the 
exact form of the discrete dispersion equation for an elec- 
tron in the adatom in the single channel case 0, 0, H, HH 
is given by 



= 0. 



(23) 



This is equivalent to a quartic dispersion polynomial after 
squaring. Note the presence of the singularities in the 
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FIG. 6: The wave numbers K± for the twelve solutions of the set of equations (|21[) . The arrows represent how the solutions 
move when we increase Ed from —2th to 2th with £ h = 0.345th and g = 0.1th- The left-hand figures (a) and (c) show the K+ 
plane and the right-hand figures (b) and (d) show the K- plane. In the top figures (a) and (b), the solutions in the upper half 
plane are in the Riemann sheet I and those in the lower half plane are in the Riemann sheet IV. In the figure (c), the solutions 
in the lower half plane are in the Riemann sheet II and those in the upper half plane are in the Riemann sheet III. In the 
figure (d), the solutions in the upper half plane are in the Riemann sheet II and those in the lower half plane are in the Riemann 
sheet=III. In the bottom figures (c) and (d), the solutions RQ4 and RQ5 cross the real axes from the Riemann sheet III into 
the Riemann sheet II when _Ed is increased from negative to positive. The vertical gray lines represent K± = — n, 0,n. 



third term at z — ±<h + t' h ; these are a result of the 
singularities in the density of states function p_ just as 
in the two-channel case. 

Now we will obtain an approximate form for the energy 
of the persistent stable state as a solution E ps of the 
single channel discrete dispersion equation (f25)) . This 
approximation will hold under the assumption that the 
impurity energy is much larger than that of the lower 
inner energy band, that is —t^ + t' h . The energy of 

this persistent stable state and its relation to the energy 
band E_ are represented diagrammatically in Fig. 1121 



the definition of E_ is given in Eq. (|7|) . Considering this 
observation, we write an expansion for E ps near the lower 
band edge — + t' h in powers of the coupling constant as 

E ps - (-< h + t' h ) + Xa g a + XW + ■ ■ ■ » ( 24 ) 

in which < a < 7, ~ 1 is independent of the cou- 
pling g at every order, and the power of g in each term is 
to be determined below. Notice that we have "skipped" 
the (3 term in Eq. (|2"4")l ; this is anticipation of the appear- 
ance of the decay rate in a similar approximation that we 
will perform for the QBIC state further on (cf. Eq. (l2"9)) ) . 
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FIG. 7: The energy for the twelve solutions of the dispersion equation, put in the corresponding Riemann sheets. The arrows 
represent how the solutions move when we increase Ed from — 2th to 2th with t' h = 0.345£h and g = 0.1th. Energy is measured 
in units of th = 1. The vertical gray lines represent the four van Hove singularities. 



We can now use Eq. (flM)) to write the dispersion equa- 
tion as 



z-E d = 



9 



1 



2 v/-2i hX «.9 Q - 2t hX ^sn + 0(5 Q+7 ,3 2Q ) 

(25) 

We can expand this to obtain 



(-t h + t' h - E d ) + Xa g a 

7 2-a/2 



(26) 



9 



2\/-2thXa 



^7-a+(2-a/2) )5 a+(2-a/2)-j_ 



The term in parentheses on the LHS and the first term 
on the RHS represent the two first-order terms; this can 
be shown to be the only consistent choice. Since the term 



in parenthesis is zeroth order in g, equating these terms 
gives the condition 2 — a/2 = 0, or a — 4, as well as 

The second-order correction is then given by equating 
the second term on the LHS of Eq. (f2"6"|) with the second 
term on the RHS; this gives 7 = 2a = 8. Making use 
of Eq. (|27[) . we can then write the expansion for the real 
energy for the persistent stable state Eq. (|24[) as 



E, 



ps 



(-*h + 4) 



8th (th ~ th — E,_\) 



7 g 4 + 0(g s ). (28) 



We see that the energy shift from the band edge —t\ l + 
t' h is of order g 4 . Notice that it was the cancellation 
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FIG. 8: Real part of the energy for the twelve solutions of the dispersion equation as a function of -Ed for the values t' h — 0.345£h 
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FIG. 9: Imaginary part of the energy for the eight complex solutions of the dispersion equation as a function of Ed for the 
values t' h = 0.345th and g — 0.1th. The unit of the energy is th = 1. The vertical gray lines represent the four van Hove 
singularities. 



of the band edge term — th + tL when we plugged the 
expansion ([24]) into Eq. ([23]) that resulted in the g 2 " Q / 2 
term in Eq. ([2"rj]k hence the persistent stable state is a 
direct result of the divergent van Hove singularity at the 
band edge. Also, note that E ps is purely real, including 
higher orders. 



2. Analytic approximation for the QBIC state Q2 in the 
two-channel model 



Now we will obtain a similar expansion for the en- 
ergy of the QBIC state Q2 in the case of the two-channel 
model that is the main subject of this paper. We again 
present a diagrammatic representation of the energy of 
this state and its relation to the two energy bands in 
Fig. [T31 As before, if we assume that 3> — th+t' h , then 



the real part of the energy of this state will be shifted such 
that it lies slightly below the lower edge of the upper en- 
ergy band at — th +t' h just as in the case of the persistent 
stable state. However, in this case the band edge (and 
therefore the energy of the state Q2 as well) is embedded 
in the continuum of the lower energy band. It is well 
known that the continuum will have a de-stabilizing ef- 
fect on a discrete state that is embedded within it. In this 
case the embedding will result in a slight de-stabilization 
of the otherwise stable state. For a further illustration of 
the relationship between the persistent stable state and 
the QBIC state, refer to Fig. fTiTa.bV 

As before, we write the expansion for the energy of the 
state Q2 as 

E Q2 = (-th + O + x a g a + xpg + x 7 .9 7 + ■ • • (29) 

where < a < f3 < 7. Here we have included the /3 
term; below we will see that this term results in the small 
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decay rate giving the QBIC effect. The zeroth-order term 
— tii + t'h places the solution well inside the continuum of 
the lower energy band. Looking to the discrete dispersion 
relation ([TO]) , we note that it is the second term in the 
square brackets that is associated with the upper energy 
band E- including the embedded band edge at — th + t' h ; 
the first term is associated with the embedding lower 
energy band E + [28| . 

Putting the expansion ([29]) into Eq. (fTO)) . we obtain 



(-t h + t' h - E d ) + Xa g° 



2-a/2 



X09 



fi-a+(2-a/2) 



+ .9' 



2-a/2 



0(g 



7— a 



9 ,9 



2(0- 



set of parentheses on the RHS. This gives the condition 
on /? as f3 — a = 2, or (3 — 6, while we also obtain 



±i 



X0 



I6t h (t' h - t h - E d ) 3 y/f h {t h -f h ) 



(31) 



This correction is purely imaginary as th — t' h > in 
the case of overlapping bands. This is the QBIC effect, 
appearing at the level of a second-order perturbation cal- 
culation. It is a direct result of the interaction of the dis- 
crete state with the two overlapping energy bands, as we 
have just shown. The remaining order terms in Eq. (|30|) 
yield the consistent result that 7 = 8 as in the case of 
the persistent stable state before. 

The expansion ([2TJ|) for the state Q2 can now be written 

as 



Eq2 — (~th + th) - 
±i 



1 



79 



(32) 



8t h (f h -t h -E d ) 2 

g 6 + 0(g 8 ). 

m h (t' h - t h - E d )3y/t h (t h -t' h ) y 



Comparison with Eq. ([28)) above emphasizes that this 
state behaves essentially like the persistent stable state 
that results from the van Hove singularity; only in this 
case the energy shift puts this state in the continuum of 
the lower energy band with a small decay rate at order 

For a numerical comparison, we can plug in the num- 
bers t' h = 0.345th, 9 = 0.1th, and E d = 0.3t h from 



Table |TJ Plugging these numbers into Eq. (|2"8")l for 
the persistent stable state in the single channel case 
gives E ps ss —0.655013701th (this value includes the g 8 
term, not explicitly given in Eq. ([2"B"|) above). Plugging 
the same values into Eq. (f3"2"|) gives the value Eq 2 rs 
-0.655013704th - £(1.5095 x lO" 7 )^ in agreement with 
the numerically obtained value reported in Table |TJ 

Finally, note that the expression for Eq 2 given above 
diverges in the case tL — th- This is an indication that 
Eq. (|32| breaks down as t' h approaches th- We will find 
a new expression to replace Eq. ([3^]) in the special case 
t' = t h in Sec. IV. 



'h 



0(g a " 



(30) 



Again we note that in this expression the terms in the 
first set of parentheses on the RHS are associated with 
the embedded singularity at — t b +t' h and those in the sec- 
ond set of parentheses are associated with the embedding 
energy band. We will see that equating one term from 
each results directly in the QBIC effect. The first-order 
correction is obtained exactly as in the case of the persis- 
tent stable state above; by equating the term — th+t' h — E d 
on the LHS with the first term in the first set of paren- 
theses we obtain the condition a = 4 and Eq. ([27]) as 
before. The second-order condition is then obtained by 
equating the second term in the first set of parentheses 
on the RHS of Eq. ([3171) with the first term in the second 



Wave function analysis and numerical 
simulations of time evolution 



In this subsection we will look more closely at the wave 
function for the QBIC state Q2. In particular, we will 
show that the wave function for Q2 appears to be lo- 
calized near x = 0, although it actually behaves as a 
decaying state with an exponential divergence for large 
x. We will also verify that the "— " channel provides 
the dominant contribution to this wave function in the 
vicinity of the origin, as would be expected since it is 
the singularity associated with this channel that results 
in the nearly localized behavior of this state. (Note that 
we have added quotation marks on the minus sign above 
to avoid any confusion in the notation. Hereafter, we 
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FIG. 11: Imaginary part of the energy (decay rate) for solution Q2 (Sheet II) as a function of .Ed for the values t' h = 0.345th 
and g = O.lih. (a) A linear plot on the and (b) a semi-logarithmic plot. The unit of the energy is th ~ 1. The vertical gray 
lines represent the four van Hove singularities, just as in Fig. [5] For the solution Q2, the decay rate is amplified (such that 
ImE ~ <? 4//3 ) in the vicinity of the lowest outer band edge — th — t' h . For values of S> th — t' h , however, the decay rate 
becomes small (with ImE ~ g G ) as the solution behaves as a QBIC state. 




FIG. 12: Diagrammatic representation of the upper energy 
band E- associated with the (shifted) single channel Hamil- 
tonian and the purely real energy associated with the 
persistent stable state E pB . 
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FIG. 13: Diagrammatic representation of the two energy 
bands E± associated with the full two-channel Hamiltonian 
Tt of Eq. JlJ and the real part of the energy associated with 
the quasi-bound state in continuum Eq2 that lies below the 
lower inner band edge at — th + t' h . 
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FIG. 14: (a) A schematic view of the strongly bound state 
(due to the van Hove singularity) of a one-channel system 
with the eigenvalue just below the lower band edge, (b) Some 
of the bound particles leak into the attached channel in the 
two-channel model. 



more gradual time scale than an ordinary decaying state. 



Wave function for QBIC state Q2 in the y — 1,2 basis 



will drop the quotation marks and simply write this as: 
— channel.) Finally, we will also compare the time evolu- 
tion of the ordinary decaying state SI with that of Q2 in 
order to demonstrate that the QBIC decays on a much 



In Fig. fTSFa.b) we plot a numerical result of the wave 
function from Eq. (116p in the non-diagonalized channels 
y = 1,2 for the state Q2. Looking at Fig. [ToT a) (linear 
scale) we see that the wave function for Q2 appears to be 
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FIG. 15: (a) The wave function modulus \ip(x, y)\ of the state 
Q2 around the origin (linear scale), (b) The same but away 
from the origin on the logarithmic scale. The (overlapping) 
plots for y — 1 (the upper leg) and y = 2 (the lower leg) are 
almost indiscernible. The parameters are set to t' h — 0.345th, 
g — 0.1th and Ed = 0.3th- The wave function is normalized 
such that ipa = 1. 



localized for values of x near the origin (where the adatom 
is attached to the wire system). However, in Fig. [T5Tb) 
(logarithmic scale) we see that the wave function indeed 
behaves as that for a decaying state with an exponential 
divergence in space far away from the origin. Also notice 
that the contributions to the wave function from the two 
original channels y = 1, 2 are almost exactly equal in 
either plot (the two graphs are overlapping). 



2. Wave function for QBIC state Q2 in the a — +, — basis 

In Fig. |T5] we plot separately a numerical result of 
the wave function contribution to the state Q2 from the 
+ and — channels (as exemplified in Eq. (fTS"]) ) in the ba- 
sis of the partially diagonalized Hamiltonian ([4} . We also 
plot the discrete wave function associated with the impu- 
rity state, each at initial time t = with the usual choice 
of parameters t' h = 0.345th, g = 0.1th and = 0.3th- 
Keeping in mind the numerical choice for the coupling 
constant g = 0.1th, we can see in this figure that at 
the origin the amplitude of the wave function |t/>_(:e)| 
for state Q2 in the — channel is of the order g larger 
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FIG. 16: The eigenfunction of the state Q2 for t h = 0.345t h , 
g = 0.1th and Ed = 0.3th on a logarithmic scale as a func- 
tion of position. The amplitude modulus of the — channel, 
|-i/)_(a;)|, that of the + channel, and that of the dot, 

\ipd\, are indicated. The wave function is normalized such that 
V>d = 1- 



than that of the adatom wave function \tpd\, and that 
in turn |^d| is order g larger than the amplitude of the 
wave function in the + channel. This implies 

that |^ + (x)| is order g 2 smaller than owing to 

the fact that the singularity in the — channel gives the 
localized behavior of the state near the origin. 

We may analytically demonstrate the relative orders of 
the wave functions by making use of the results that we 
obtained for the QBIC state previously in Sec. II-B. First 
we may add or subtract the first two equations in (fT7|) to 
obtain 



l^dl = \/2th, 

^±(0)1 g 1 

and the channel weight function 



sin K ± | 



(33) 



0(E) 



\M0)\ 
^-(0)| 



Isin^l _ Jt h 2 -(E-t h )* 



! sin K< 



th 2 -(E + t' h y 



(34) 



in which we have made use of Eq. (fl2|) to write ^±(0) = 
A± at the origin x = 0. These relations hold for any of 
the twelve eigenstates of the ladder system. 

We here use the relations for the QBIC eigenvalue 
E = Eq2 and the corresponding wave numbers K±. We 
can now use the expansion (|32[) with the continuous dis- 
persion equations ([7]) to obtain 



*jr + - a J4(i-4 ] h 



and 



sini^_ = i- 



2th (t'h — th — Ed 
Applying these in Eq. (f3"3"| gives 

IV> d | 



■9 2 + 0(g 4 ). 



I^±(0)| 



(35) 



(36) 



(37) 
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in agreement with our discussion of Fig. [16] above. Fi- 
nally, the channel weight function (|34[) for the state Q2 
is given by 

e(E ) ^ +(0)l = gHh2 

Q l^-(0)| 4 v ^(V r t[)|-th + 4- E d \' 

(38) 

such that the contribution of the + channel is order g 2 
smaller than that of the — channel. 



3. Generic channel weight function 

We may generalize the preceding discussion to states 
other than Q2. As an application of Eq. ([3"4"| . consider 
the decaying states SI and Q2 in the vicinity of the outer 
band edge at — ij, — t' h . Here the state Q2 will not behave 
as a QBIC state with a small decay rate, but instead has 
an amplified decay rate of order g 4 ^ 3 in the vicinity of 
the singularity at the outer band edge (for example, see 
Fig. [TIT a)) J2j. The state SI also has an amplified decay 
rate for this range Ed ~ — — t' h . For either of these so- 
lutions, Eq. ([3"4"[) gives 0(£'si.Q2) ~ 9~ 2 ^ 3 , demonstrating 
that in this case, it is the lower channel + that provides 
the largest contribution to the wave function. The reason 
for this is that it is the van Hove singularity at the outer 
band edge —th — t' h (associated with the lower band edge 
E + ) that results in the amplification of the decay rate, 
while the upper channel E_ plays little role in this effect 
in this region of the energy spectrum. 

4- Time evolution for QBIC state Q2 against ordinary 
decaying state SI 

In Fig. [T7] we show the time evolution &(x, y, t) and 
^d(i) for the states SI and Q2 for the ordinary choice of 
the parameters t' h — 0.345ih, 9 — O.lfh an( l £?d = 0.3th, 
under which the state Q2 will behave as a QBIC state. 
We see that indeed on the time scale under which the 
ordinary state SI decays almost completely, the state Q2 
appears to behave as a localized state without a notice- 
able decay rate (for this time scale). The details of the 
numerical method by which we have obtained the time 
evolution simulations in these plots are presented in Ap- 
pendix A. 



case t' h — th- Here we will find that the 12th-order dis- 
persion polynomial reduces to a lOth-order polynomial, 
while the QBIC decay rate is amplified such that it is 
proportional to g to first order; both embedded singu- 
larities play a role in this modified QBIC effect. 



A. Energy spectrum for t' h — case 

Here we will comment on the case tl = 0, in which the 
chain-to-chain hopping parameter vanishes. Since the 
adatom is coupled to only one chain, it is to be expected 
that this system should reduce to the single chain model 
(along with an additional uncoupled chain). Indeed, con- 
sidering the discrete dispersion equation ([TO)) we see that 
for t' h — the two terms on the RHS containing the 
square roots will agree. Therefore, in the case where the 
sign of the two square roots agrees, these two terms will 
combine to give a new dispersion equation equivalent to 
that for the single-channel model (equal to Eq. ([23)) after 
setting t' h — and replacing with the re-normalized cou- 
pling constant g' —> \[2g ) . The case in which the sign of 
these two terms agree corresponds to Sheets I and IV in 
the complex energy surface, which contain four solutions 
(two in each sheet). These four solutions then behave 
precisely as the four solutions to the quartic dispersion 
polynomial in the original single-channel model. 

Meanwhile, for the case in which the sign of the two 
square roots in Eq. (flU)) are opposite, then these two 
terms cancel and the dispersion equation becomes trivial. 
This case corresponds to Sheets II and III, which contain 
eight solutions. Meanwhile, since the two branch cuts in 
this case overlap exactly, Sheets II and III become math- 
ematically (and physically) inaccessible. If one travels 
through the branch cut in Sheet I, one will always appear 
in Sheet IV (and vice versa). Hence the energy surface is 
effectively two-sheeted and essentially equivalent to that 
in the single chain system. 



B. Energy spectrum for t'h — th case 



IV. ENERGY SPECTRUM ANALYSIS IN TWO 
SPECIAL CASES 

In this section, we will briefly consider the energy spec- 
trum analysis for two special cases. In the first special 
case (t' h = 0) the two-channel model reduces to the single 
channel model when the chain-to-chain hopping param- 
eter t' h vanishes. Even though this case is trivial, it is 
instructive to see the relation of the two-channel model 
to the single channel model. Then we will consider the 



In the case t' h — t^, the two inner band edges will 
overlap to form a single embedded band edge as indicated 
in Fig.[l8j Note that, in a certain sense, both band edges 
are still present. In fact, we will find that in this case 
the van Hove singularity from one overlapping band edge 
will actually amplify the QBIC decay rate that results 
from the van Hove singularity of the other band edge, 
so that the decay rate will now be proportional to g 4 . 
This is a unique combination of two previous effects, both 
resulting from the van Hove singularity. 
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Solution Q2 for E d = 0.3t h 



FIG. 17: Time dependence of wave function of solutions SI ((a)-(d)) and Q2 ((e)-(h)) for t' h — 0.345th, g = 0.1th and 
Ed = 0.3th- The solid (red) curves represent ^(a;, l,t), the broken (green) curves represent ^(x,2, t) and the dots represent 
^d(t)- The wave functions are normalized such that ^d(O) = 1. The QBIC solution Q2 clearly decays on a much slower time 
scale than that of the ordinary decay state SI. 
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to the general case. Specifically, the two solutions RQ4 
and RQ5 are no longer present. Note that the QBIC 
solutions Q2 and R2 remain in this simplified system. 

As was done in Sec. Ill A and in Fig. [51 it is more 
convenient for the purpose of placement of each solution 
in the correct Riemann sheet to solve 



FIG. 18: Band structure for the t h = th case with a single 
"unified" embedded band edge in the center of the spectrum. 



1. Tenth order dispersion polynomial for the t' h — th case 

If we set t' h — ^ in Eq. (flO|) , then the discrete disper- 
sion equation becomes 



E d 



1 



1 



+ 2t h y/z - 2t h 



0. (39) 



Note that here the position of a solution in the com- 
plex energy surface should be determined only by the 
two square roots 1/y/z ± 2th! the overall factor of 
plays no role in this determination. However, this overall 
factor represents the presence of a van Hove singularity 
at z = in both bands. 

As before, we can square this equation twice to find an 
equivalent tenth-order dispersion polynomial; hence two 
solutions have vanished from the system in comparison 



-th cos K + — th = — ih cos K- + th 
1 1 



= E d 



9 



2ith sin K + 2ith sin K_ 



(40) 



with respect to K± than to solve the discrete dispersion 
equation ()39() with respect to z — E. The imaginary 
parts of the solutions of the above simultaneous equations 
give the correct Riemann sheet. The eigenenergy of each 
solution is given by the continuous dispersion equation 
E = — th cos K± =F th- 

The dependence of the real and imaginary parts of all 
the solutions on E d can be seen in Figs. HU and [20l re- 
spectively, for the value g = 0.1th for the coupling. 



Modified QBIC effect 



As in the general case, the solutions Q2 (for E& > th) 
and R2 (for E& < — th) arc QBIC states. Their real parts 
behave similar to the persistent stable states in that they 
lie near to the embedded band edge at z = in Fig. [TIB 
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FIG. 19: Real part of the energy for the ten solutions of the simplified dispersion equation (|39[1 as a function of for the 
special case t' h — th with the choice g = 0.1th- The unit of the energy is — 1. In the top four figures (a)-(d), each solution 
is plotted in the corresponding Riemann sheet. The overlapping curves represent a complex conjugate pair, for which the real 
part of the energy is exactly the same. The vertical and horizontal gray lines represent the three van Hove singularities. In the 
bottom figure (e), a part of the top four figures (a)-(d) are shown with all Riemann sheets superimposed. 
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FIG. 20: Imaginary part of the energy for the six complex solutions of the simplified dispersion equation (|39[1 as a function of 
Ed for the special case t'^ = ih with the choice g = 0.1th- The unit of the energy is th = 1. The vertical gray lines represent the 
three van Hove singularities. 



However, they have a small non-zero decay rate as can 
be seen in Fig. [201 We can obtain an approximate form 
for the energy eigenvalue for these states in a manner 
similar to that we employed before in the general case. 
Figure fT9l shows that the real part of the energy for Q2 is 
also small around the origin. Hence, we assume that the 
expansion of Eq2 begins with the order g a with a > as 
Eq2 = Xa9 a + • • • • Applying this expansion in Eq. (|59"|) 
and using a similar argument as that given in Sec. Ill B 
gives 

fxZ=9 2 - a ' \ (41) 



2y / 2t} 1 Ea 

in which the two bracketed terms (both of which are nec- 
essary for Eq2 to be imaginary) have been contributed re- 
spectively by the two bracketed terms on the LHS of (f59")) . 
and the factor 1 / y/\z\ associated with the van Hove sin- 
gularity in both terms has resulted in the factor g 2 ~ a l 2 . 



We then obtain the familiar condition for the order of a 
as a = 4 as well as the condition \a — —i/i^hE^), from 
which we obtain 



EQ2 — 



it h E 2 



(42) 



We see that the state Q2 is still quasi-stable in compari- 
son to the ordinary decay rate proportional to g 2 . How- 
ever, the QBIC decay rate has been amplified, so that it 
is proportional to g 4 (instead of g e in the general case) 
due to the overlapping band edges. In a sense, the em- 
bedded van Hove singularity at z = from the £L band 
has resulted in a QBIC state while the singularity from 
the E + band has amplified the usual decay rate of g e to 
,9 4 - 



20 



Indeed, the channel weight function (|34j) reduces to 

«^ = ^ <«> 

for tL = th and gives 6{Eq?) ~ i for the specific case of 
the QBIC state. Thus both channels contribute equally 
in this modified QBIC effect. 

V. CONCLUDING REMARKS 

Before we conclude this paper let us give a brief com- 
ment on the order of the dispersion polynomial and num- 
ber of solutions (and QBIC solutions) for an n-channel 
wire model. In the single channel model 0, @, [13] the 
quartic dispersion polynomial yielded four solutions. If 
we consider the position of the real part of these solu- 
tions (see, for example, Fig. 2(a) in [2|), then regardless 
of the value of Ed there are always two solutions that lie 
near the value Ed and one solution that lies near each of 
the two band edges; when Ed lies well within the con- 
duction band, both of the latter solutions are persistent 
stable states that co-exist with the decay solution. We 
can conceptualize this system with the statement that we 
have 2e a + 2be = 4 solutions, where 2^ d indicates that 
there are two solutions near Ed and 2be indicates that 
there are two more solutions near the two band edges of 
a channel. 

In general, there are two band edges associated with 
each channel and hence the total number of the band 
edges is TVbe = 2n. In the present case we can see in 
Fig. [8] that for the two-channel model there are always 
four solutions that lie near the energy Ed {^E d ), while 
there are two solutions near each of the four band edges 
(2x4 B e); hence we have 4 Bd +2x4 B E = 12 solutions. We 
can generalize this observation by saying that for an n- 
channel model there are (2 n ).g d solutions associated with 
the discrete energy Ed and 2™ _1 x 2Vbe = 2"n solutions 
associated with the -/Vbe band edges. Hence there are 
-'Vsoins = 2"(1 + n) total solutions for an n-channel quan- 
tum wire coupled with a single discrete state, given as 
the solutions to an A^ so i ns -order dispersion polynomial. 
As mentioned in Sec. II B, these solutions will live in a 
Riemann surface composed of 2n Riemann sheets. 

Regarding the number of QBIC states for the n- 
channcl model, this number will depend not only on the 
number of embedded band edges, but also on the value 
of Ed- However, we can state that in general the maxi- 
mum number of QBIC states for each model is equal to 
the number of band edges which could possibly be em- 
bedded in the continuum of another band. This number 
is equal to 2 n_1 , although this equation is obviously not 
valid in the single-channel model for which no embed- 
ding is possible. We intend to present these statements 
in greater detail elsewhere. 

We have demonstrated the existence of the QBIC state 
in the context of a two-channel quantum wire with an at- 
tached adatom impurity. Ordinarily, an electron in the 



impurity with an energy deep inside of the conduction 
band of the wire would be expected to decay and travel 
along the length of the wire. We would also expect that 
we could describe the decay rate using Fermi's golden 
rule. However, due to the combined effect of two over- 
lapping conduction bands (with van Hove singularities at 
the band edges) , we have shown that a QBIC electron will 
remain stable inside the impurity for ordinary time scales. 
In particular, we have demonstrated the connection be- 
tween the QBIC state and the persistent stable state that 
results from the van Hove singularity at the band edge 
in the single channel model. In the two-channel model, 
this persistent stable state is slightly de-stabilized by the 
presence of a second conduction band. 

We have also shown in Eq. (|52")l that the characteristic 
decay rate for the QBIC is generally on the order of g 6 , 
although this may be modified under certain conditions, 
such as the case t' h = th in the two-channel model (under 
which ImE ~ g ). 

While we have demonstrated the above specifically for 
the two-channel quantum wire, it is easy to show that 
this effect should occur in other one-dimensional models 
which have the characteristic square root divergence in 
the DOS given in Eq. ([8]). This includes models such as 
an electromagnetic waveguide [5| when we consider two 
overlapping field modes, as we mentioned above. Hence, 
the origin of the QBIC effect is quite different than that 
of the BIC effect originally proposed by von Neumann 
and Wigner. While we can associate each QBIC state 
with a divergent band edge singularity embedded in the 
continuum of another energy band, the BIC states are 
associated with zeros in the interaction potential that oc- 
cur in the continuous energy spectrum for certain models 
with an oscillating potential. Hence the QBIC is more 
closely associated with the DOS function while the BIC 
is more closely associated with the interaction poten- 
tial. It is also conceivable that the QBIC effect may 
appear in some two-dimensional systems (such as a two- 
dimensional tight-binding lattice) that have a character- 
istic logarithmic divergence in the DOS 

As we remarked in our previous Letter [l[ , because our 
quasibound state has a small decay rate (imaginary com- 
ponent of the eigenenergy) , it is not, strictly speaking, "in 
continuum." However, this decay rate is extremely small, 
such that the QBIC state should behave as if it were a 
bound state with real part of the eigenenergy deeply em- 
bedded in the continuum even on relatively large time 
scales. In this sense, the QBIC will essentially behave 
as the BIC under actual experimental conditions. Mean- 
while, the BIC is a true bound state with a purely real 
energy spectrum, but only under ideal conditions. Since 
the BIC exists only at discrete points (with zero mea- 
sure) in parameter space any noise in the system (such 
as thermal noise) in an experiment may actually lead to 
a small decay rate for the BIC. On the contrary, since the 
QBIC exists for a wide range of parameter space (with 
non-zero measure), it is robust against noise. It may be 
easier to prepare the QBIC in the experiment. 
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It should be mentioned that there may be certain mod- 
els with embedded singularities in which the QBIC effect 
due to the DOS singularities will be washed out as a re- 
sult of the form of the interaction potential. There is 
at least one example of a single channel model in which 
the interaction potential washes out the effects of the 
singularity and prevents the persistent stable state from 
forming [181 ]. 
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where y — 1,2. The vector form of the wave functions is 
given by 



e- iB %(x). (A4) 



We restrict the region \x\ < L, (L > 1) to compute the 
time evolution of the wave function. The steady wave 
function ip( x ) m Eq. (| 16[) has the following recursion 
property: 



$ y (L+l) = U$ a (L + l) 

= u( elK Q + f K _ \ u- l u ML) 

= V cS $ V {L), (A5) 



where we have used Eqs. (|13l) and (I16p . We have defined 
the effective potential V e s as 
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APPENDIX A: NUMERICAL METHOD FOR 
TIME EVOLUTION SIMULATION OF WAVE 
FUNCTIONS 

In this appendix we describe our numerical method for 
obtaining the time evolution of the wave functions of the 
resonant states of the two-channel Hamiltonian. We rely 
on the method which was proposed in previous work [221 ] 
to solve the time-dependent Schrodinger equation accu- 
rately (despite truncation of the domain of x in numerical 
calculations). The time-dependent Schrodinger equation 
is given by 

ift^|*(t)> =#|*(t)>, (Al) 
and the initial condition is fixed as 

|*(0)) = (A2) 

where \^) is the eigenstate of H in Eq. (flT]) . We define 
the time-dependent wave function as follows: 



which we have extended from the scalar form of the chain 
model @, [H, [H[. At x = ±L (L > 0), using the effec- 
tive potential (|A6|) , the left term of the time-dependent 
Schrodinger equation (|A1[) becomes 

= -| {$ y (±(L - 1), t) + $ y (±(L + 1), t)} 

= -|^(±(L-l),<)-|v/ cff ^(±L,i), 

-t' h S%(±L,t), (A7) 

in which we define the matrix S by 



* d (t) = <d|*(t)), *(x,y;t) = (x,y\*(t)), (A3) 



Thus, we can obtain the time-dependent Schrodinger 
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equation in the closed region \x\ < L as follows: and 

'- t -±-{$ y (±(L-l),t) + V eS $ y {±L,t)} 
-t' h S%(±L,t) 

for x = ±L, 

-%{$ y (x-l,t)+$ y (x + l,t)} ih^-y d (t)=E d y d (t)+gV(0,l;t). (A10) 

-t' h S$ y (x,t) 

for 1 < \x\ < L - 1, 

-^{# y (-M) + # y (i,t)} 

-t' h S %(0,t)+g* d (t) Q 

for x — 0, This is the method by which we have produced the time 

(A9) evolution simulations presented in Fig. [T7l 
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